Existence of hysteresis in the Kuramoto model with bimodal frequency distributions 
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We investigate the transition to synchronization in the Kuramoto model with bimodal distributions of the 
natural frequencies. Previous studies have concluded that the model exhibits a hysteretic phase transition if 
the bimodal distribution is close to a unimodal one, due to the shallowness the central dip. Here we show that 
proximity to the unimodal-bimodal border does not necessarily imply hysteresis when the width, but not the 
depth, of the central dip tends to zero. We draw this conclusion from a detailed study of the Kuramoto model 
with a suitable family of bimodal distributions. 

PACS numbers: 05.45.Xt 



I. INTRODUCTION 

Understanding the dynamics of large populations of het- 
erogeneous self-sustained oscillatory units is of great inter- 
est because they occur in a wide range of natural phenomena 
and technological applications |1]. Often a macroscopic sys- 
tem self-organizes into a synchronous state, in which a certain 
fraction of its units acquires a common frequency. This occurs 
as a consequence of the mutual interactions among the oscil- 
lators and despite the differences in their rhythms [2]. Exam- 
ples of collective synchronization include pacemaker cells in 
the heart and nervous system 01 > synchronously flashing 
fireflies [5], collective oscillations of pancreatic beta cells |^ 
and pedestrian induced oscillations in bridges [7]. 

A fundamental contribution to the study of collective syn- 
chronization was the model proposed by Kuramoto [St]. This 
model, and a large number of extensions of it, has been ex- 
tensively studied because it is analytically tractable but still 
captures the essential dynamics of collective synchronization 
phenomena (for reviews see fll|9l[l3,llitl)- The original Ku- 
ramoto model consists of a population of N oscillators inter- 
acting all to all. The state of an oscillator i is described by its 
phase di {t) that evolves in time according to 



K ^ 



(1) 



The parameter K determines the strength of the interaction 
between one oscillator and another. The oscillators are con- 
sidered to have different natural frequencies oji, that are taken 
from a probability distribution g{uj). In his analysis Ku- 
ramoto adopted the thermodynamic limit N ^ oo and con- 
sidered g{u!) to be symmetric. In this case, and without 
loss of generality, the distribution can always be centered at 
zero, i.e. g{uj) — g{—uj), by going into a rotating framework 
Oj — >■ 9j -\- ^t. 

Kuramoto found useful to study the synchronization dy- 
namics of system ([TJ in terms of a complex order parame- 
ter z = N^^ ^f=i exp(z 6j). Note that z is a mean field 
that indicates the onset of coherence due to synchronization in 



the population. System ((TJ possesses an incoherent state with 
z — (that exists for all values of the coupling strength K) 
in which the oscillators rotate independently as if they were 
uncoupled, di (t) ^ Uit. Using a self-consistency argument, 
Kuramoto found that for a unimodal distribution g{Lo), above 
the coupling's critical value 

Kc^Arry (2) 

7r.g(0) 

a new solution with asymptotics 



' K-K, 
-7rg"(0) 



(3) 



branches off the incoherent (z = 0) solution. This emerg- 
ing solution is a partially synchronized (PS) state, in which 
a subset of the population S entrains to the central frequency 
{9ies = const.). 

Equation (O shows that the orientation of the PS bifurcat- 
ing branch depends on whether the distribution is concave or 
convex at its center As a consequence of that, at K = Kc 
the PS state is expected to bifurcate supercritically for uni- 
modal distributions (.g"(0) < 0) and subcritically for bimodal 
distributions (.g"(0) > 0). However, Kuramoto's analysis did 
not permit to study the stability of the solutions and thus one 
cannot conclude whether bimodal distributions show bistabil- 
ity close to the transition point ^ (see discussion in p. 75 of 
Jstl). In fact, Kuramoto discarded the possibility of bistabil- 
ity. Instead he expected the incoherent state to become un- 
stable earlier, i.e. at a certain critical value K'^ < K^, via 
the formation of two symmetric clusters of synchronized os- 
cillators near the distribution's maxima (later Crawford called 
this state standing wave (SW) |12]). As the coupling is in- 
creased further, he predicted that the interaction between the 
clusters would tend to synchronize them forming a single syn- 
chronized group, i.e. a PS state. 



A. Sum of unimodal distributions with different mean 

After Kuramoto's seminal work [8], several articles have 
further investigated the synchronization transition in model 
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FIG. 1 : (Color online) Examples of bimodal frequency distributions 
given by Eq. ^ with 5=1. Left panel: ^ = 7 (what implies 
g{0) = 0). Note that as 7 decreases the maxima of the distribu- 
tion become closer For all these distributions (with C — 7) the 
route to synchronization as K is increased from zero is I— >SW-^PS, 
c.f. Fig. [3] Right panel: Two examples with ^ < 7. The distribution 
depicted with a continuous line has well separated peaks and shows 
a transition I^SW^PS, whereas the other distribution is closer to 
the unimodal limit ^ and presents hysteresis in the route to synchro- 
nization, c.f. Fig.|4] 



dIJ with symmetric bimodal distributions IH, [Tsi [TillTsi [Tel 
[l7il . These studies assumed g{Lu) to be the superposition of 
two identical even unimodal distributions (i{uj) centered at 
±0^0: g{i^) ~ g{oj + ijJq) + g{Lu — ujq) |32]. Parameter o^o 
controls the separation of the peaks. Decreasing ojq the distri- 
bution's maxima approach each other and, at the same time, 
the central distribution's dip becomes shallower (i.e. g{0) in- 
creases). Eventually, at a value lu^ — luqb that satisfies 



= 0. 



(4) 



the peaks merge and the distribution becomes unimodal. The 
dynamics of the Kuramoto model for distributions of this type 
is as follows 1 13]: When the peaks are well separated (ujq 
larger than a certain value luqd) the transitions increasing K 
are as Kuramoto foresaw: Incoherence SW — > PS. How- 
ever, if the peaks are near (lucid > ^^q > ^ob) there exists a 
range of K below Kc where bistability between incoherence 
and either a PS or a SW state is observed, as Eq. (O suggested 



B. Difference of unimodal distributions with different width 

In this article we are interested in understanding the syn- 
chronization transition in the Kuramoto model with bimodal 
distributions in situations that cannot be achieved summing 
even unimodal distributions. In particular summing even dis- 
tributions implies that if the peaks are brought closer the cen- 
tral dip becomes less deep (unless the distributions are Dirac 
deltas). Thus we cannot approach the peaks arbitrarily near 
while keeping the central dip's depth (see e.g. in the left panel 
of Fig. [T] for a distribution family with constant depth but ar- 
bitrary distance between the peaks). 



We will use a family of bimodal distributions that are con- 
structed as the difference of two unimodal even functions with 
the same mean and different widths: g{ijj) ~ gi{o-!) — 52 (t^)- 
These distributions could be useful to model systems in which 
a fraction of the central natural frequencies of a population gi 
is missing due to for example, some resonance, symmetry, or 
external disturbance. 

We choose the functions gi to be Lorentzians, because of 
their mathematical tractability. Assuming 6 > j the normal- 
ized distribution reads 



5^ 



7 



(5) 



with ^ < 7 to be well defined, and S = l/((5 — ^) is the 
normalization constant. Without loss of generality we assume 
8 — \ hereafter, because this can be always achieved rescaling 
uj, time and the parameters: lo' — u)/6,t' — tS, K' — K/5, 
7' — 7/(5 and ^' = ^/(5. We will also drop the primes to 
lighten the notation. Figure [T] shows several examples of dis- 
tributions (|5]l. Distribution family (|5]l can exhibit an arbitrarily 
deep minimum while keeping the maxima as near as wished. 

The left panel of Fig. [T] shows two examples for the case 
^ = 7, which will be analyzed in detail below. This case 
implies g(0) = 0, which corresponds to the maximal value 
of the ratio ^/7 = 1. As 7 0, the central dip becomes 
infinitely naiTow and at 7 = the distribution becomes uni- 
modal. This unimodal transition is therefore discontinuous 
and satisfies 1.34.1: 



lim g" (lo = 0) = 00, 



(6) 



In addition, distribution Q also presents the regular 
unimodal-bimodal border via .g"(0) = at 



=7' 

with 7 7^ (line B in Fig.|2]l. 



(7) 



The outline of the paper is as follows: Section |ll] summa- 
rizes recent theoretical results that permit to reduce the Ku- 
ramoto model to a system of ordinary differential equations 
with complex variables. These results are then used to find 
the two ODEs that describe the dynamics of the Kuramoto 
model with distribution (|5]l. In Sec. [Ill] we study the spe- 
cial case ^ = 7, and we show that there indeed exists a 
transition to synchronization in absence of hysteresis inde- 
pendent on the separation between the distribution's maxima. 
Namely, in this case the route to synchronization is always: 
Incoherence— >SW^PS. In Sec.|IV]we study the most general 
case g(0) > 0, and determine the disposition of the differ- 
ent synchronization scenarios with respect to the unimodal- 
bimodal border. 



II. LOW DIMENSIONAL DESCRIPTION OF THE 
KURAMOTO MODEL 

We start considering the thermodynamic limit A'^ ^ 00 of 
model ([T]i. We drop hence the indices in Eq. ([T]i and introduce 
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FIG. 2: The parameter space of distribution l|5j [not defined above 
thie bisectrix ^ = 7 neither at point (1,1)]. Function lO is unimodal 
below line B and bimodal above it (shaded regions). Three lines sig- 
nal the loci of codimension-two bifurcations (A, B, and D) projected 
on the (7, ^) plane. Between lines D and B (dark grey region) the 
transition to synchronization involves hysteresis. 

the probability density for the phases /(6', i) J^HU]. Then 
/(6', uj, t) d9 duo represents the ratio of oscillators with phases 
between 9 and 6 + dO, and natural frequencies between w and 
ijj + du). The density function / obeys the continuity equation 

df _ djfv) 
dt 



(8) 



de ' 

where, the angular velocity of the oscillators v is given by 

v{9,uj,t) = uj - K f {9', Lu, t) sin{e - 9')d0' (9) 
Jo 

In the continuous formalism, the complex order parameter de- 
fined by Kuramoto becomes 



z{t) 



e'" f{9,Lo,t) d9dLo. 



(10) 



Since the density function f{9, uj, t) is real and 27r periodic in 
the 9 variable, it admits the Fourier expansion 



2n 



l + ^(/«(^,i)e™'+c.c.; 



, (11) 



where /„ — Note that the order parameter (fTOl l now 

reduces to 



g{uj)fi{uj,t)duj. 



(12) 



Substituting the Fourier series (fTTT i into the continuity equa- 
tion (ISj, and using Eq. (fT2t one gets an infinite set of integro- 
differential equations for the Fourier modes 

nK 

f„ = -inujfn + — iz*.fn-l - Zfn+l) ■ (13) 



Rec ently Ott and Antonsen (OA) found a very remarkable re- 
sult 1I22II : The ansatz 



(14) 



is a particular -and usually the asymptotic- solution of the 
infinite set of Eqs. ( fT3] l if a satisfies 



a = —luja + (-^ 



za 



(15) 



Equation (flSl l reduces to a finite set of ODEs for distributions 
g{uj) with a finite set of simple poles out of the real axis. 
Recalling fi — a the order parameter can be calculated by 
extending the integral in (fT2l i to a contour integration in the 
complex plane. This is possible since a has an analytic con- 
tinuation in the lower half a;-plane |22]. In turn only the val- 
ues of a at the poles of g{'jj) with negative imaginary part are 
relevant. 

Several recent studies show that the ansatz ( fT4l l yields pre- 
dictions in agreement with numerical simulations 113112211231 
m, HI In addition Ott and Antonsen theoretically 

support the validity of their ansatz for the case of a Lorentzian 
distribution [ 28] . So far, disagreement between the OA ansatz 
and numerical results has been shown for frequency distribu- 
tions with no spread and non-odd-symmetric coupling func- 
tion. This entails the freedom to select arbitrary values for 
some constants of motion 12911. 



A. Main Equations 

In this section we use the OA ansatz considering frequency 
distribution (|5]). This yields two ODEs governing the dynam- 
ics inside the low-dimensional OA manifold. First of all, it is 
convenient to express (HJ in partial fractions: 



1 



1 



2711 — i uj + i LU — '-fi oj + 'fi 
Then, according to Eq. (fT2] | the order parameter reads 

z*it)^E[aiit)~Ca2{t)], 



(16) 



(17) 



with ai{t) — a{uj = —i,t), and a2{t) — a{uj — ~ij,t). 
Using ( fTTb in Eq. ( fTSl ), we obtain the following two ODEs 
with complex variables that govern the evolution of the order 
parameter ( fTTj i 



di = —ai + k{ai — ^02) — k{al — ^0:2)0;^ 



012 



— —ja2 + k{ai — ^Q!2) — k{al — £,a2)al, 



(18a) 



(18b) 



with k = 'E.K/2. The phase space of Eqs. ( fTSl l is four 
dimensional, but due to the global phase shift invariance 
{ai, 012) — > (aie*^, 0126*'^) the dynamics is actually three di- 
mensional [see also Eqs. ( lAll ) in Appendix A]. 



4 



B. Fixed points 

According to Eq. ( [TtI i. the fixed points of Eqs. ( fTSl l corre- 
spond to steady states of the order parameter z. The trivial 
solution = a2 = yields z = 0, corresponding to the 
incoherent state. 

In order to calculate the non-trivial fixed points, note first 
that invariance under the action of the global rotation e*'^ al- 
lows us to choose a\—x\+ iy\ real, i.e. ol\ — x\. It follows 
from Eq. (llSal i that the fixed points lie on the subspace where 
0,2 is real too. We can therefore take ai and 012 as real (keep- 
ing in mind that a continuous of fixed points is generated un- 
der the action the neutral rotation e'*^). Hence, the equations 
for the fixed points are: 

Q = -xx^k{xx-ix2){\-x\) (19a) 



{^^--iX2^k{xx-ix2){\-x\) (19b) 



implies |q!2| > 1, and according to Eq. ( fT4] i the Fourier series 
of the density function tj, t) is divergent at cj = —17. 

We will see below that for large enough values of K there 
exist two more real solutions of P{X): < X(2) < 

1 ~ 1/k. In this case (except when becomes negative) 
such solutions indeed correspond to PS states of the original 
Kuramoto model ([T]i. 



III. BIMODAL DISTRIBUTIONS VANISHING AT THEIR 
CENTER (5 = 7) 

In this section we consider C = 7 what implies that distri- 
bution ^ vanishes at its center, g{0) — 0. In this case 7 (or 
^) becomes the parameter controlling the width of the central 
dip of g{uj), and the maxima of the distribution are located at 
(see Fig.[T] left panel): 



±7. 



(23) 



Additionally, note that these equations are symmetric under 
the reflection (a;i, X2) {—xi, —2:2). This implies that the 
solutions (with the exception of the solution at the origin) exist 
always in pairs with opposite signs (±xi, ±3:2). 



Subtracting Eq. (|19a| l from Eq. ( |19bt multiplied by -, we 



obtain X2 = 



- — 1] . This can be substituted back 



into Eq. ( I19al i to get a cubic equation in X 



^2. 



PiX) 



- k [{2k - 1)(1 - 7^) - 1 + ka^ - 7)] X' 

+ [(A:2-2fc)(l-70 + l + 2A:2e(e-7)]^ 

- kCb + kiC-l)]=0 (20) 



Each of the solutions of this equation yields two twin solutions 
with coordinates 



Xi 



±VX ^j;2=a;i[l 



k{l-X)i 



(21) 



After some algebra we obtain the relation of the solutions with 
order pai^ameter: 



2^VX 
K{l-X) 



(22) 



A steady state (a;i, X2) results in a time-independent value of 
z and hence it should correspond to a partially synchronized 
state. However, note that X can only take values within the 

range X e [0, 1 - 2^{^J +T - -^)] to have a z value 
consistent with its definition, i.e. \z\ e [0, 1]. 

As the polynomial in Eq. (|20] | is cubic, there is one real so- 
lution, ^(3), for all the parameters values. This solution lays 
in the range [0, 1] (for fc > 1 a better bound is [1 — 1/fc, 1], 
since P(l - 1/fc) = -^^ < q ^^^^ p^^) = 1 > 0). How- 
ever, it turns out that the fixed points associated to X(3) are 
'unphysical' (even though in some parameter ranges \z\ < 1). 
The reason is that the X2 coordinate, corresponding to the so- 
lution -'^(a), is always larger than 1 in absolute value. This 



A. Stability of the incoherent state 

In the incoherent state the oscillators are uniformly dis- 
tributed in the interval [Q,2tt), and thus the order parameter 
vanishes. This state corresponds to the fixed point at the ori- 
gin ai = ^2 = 0. A Unear stability analysis of Eqs. ( fTSl ) 
reveals that this fixed point undergoes a degenerate Hopf bi- 
furcation at kn = (1 + 7) /(I — 7). In terms of the original 
coupling constant K, we find 



Kh = 2 + 27. 



(24) 



A 



3.4 ^ 
0, the 



At this point the eigenvalues are imaginary Ai 2 
iy/^ and two-fold degenerate. Observe that as 7 
critical coupling for a (unimodal) Lorentzian distribution of 
unit width is recovered: Kh{i ^ 0) = Kc = 2/{'Kg{0)) = 
2. Figure [3] shows the boundary Kh in the (7, K) plane. As 
expected, we find that as the central dip of the distribution 
broadens (increasing 7) the stability region of the incoherent 
state grows. 



B. Saddle-node bifurcation 

The cubic equation (|20] i for the non-trivial fixed points be- 
comes greately simplified under the assumption C 7- 

- 72)^3 - k [{2k - 1)(1 - 7^) - 1] X^ 

[(fc2-2fc)(l-72) + l]X-72fc = 0. (25) 



Q{X) 



For 7 = the central dip vanishes, and we recover the so- 
lutions for a Lorentzian distribution X = Q,l — 1/k. When 
7 > there is a saddle-node bifurcation at fc = ksN, i-e. there 
is a transition from one (for k < ksn) to three solutions (for 
k > ksN)- ksN and 7 can be related imposing the condi- 
tion that the discriminant of Q{X) vanishes. This gives the 
following relation: 



7 



' 4 



(1 



20fc|^ 



1 



SksN{ksN + ly 



(26) 
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FIG. 3: Phase diagram for ^ = 7. For this case the synchroniza- 
tion transition never involves hysteresis. The solid lines mark the 
saddle-node (SNIC) [from Eq. ^] and the Hopf [Eq. ^] bifur- 
cations. Symbols correspond to the numerical estimation of the bi- 
furcation lines via numerical integration of the original Eq. (O with 
iV = 2000. 

There are two important asymptotic values for this bifurca- 
tion line, which expressed in terms of the original coupling 
constant K are 

KsNil^0)^2 + 6[iy^\oij), (27) 
i^SAr(7^1)-(3 + V8)(l-i^). (28) 

When K increases above Ksn the born solutions depart from 
each other X(^2) — ^{i) ^ VK — Ksn + h.o.t. One so- 
lution becomes progressively smaller (dX(i)(i^)/dX < 0), 
whereas the second one grows {dX (2){K) / dK > 0). The 
latter solution X(^2) yields a monotonically growing value of 
\z\ with K. This is not surprising because in the Kuramoto 
model, at large values of K, there exists always a stable PS 
solution with d\z\/dK > (and limx^oo \z\ = 1, i.e. full 
synchronization). We advance that the corresponding twin 
fixed points from X(^2) ^re stable, whereas the fixed points 
corresponding to are saddle. 

C. Numerical simulations and phase diagram 

In this section we construct the phase diagram with the loci 
of Hopf and saddle-node bifurcations that we have obtained 
above. Numerical simulations of the reduced Eqs. (fTSl ) were 
carried out and compared with the full model This permits 
to relate the dynamics of the variables ai 2 with the actual 
dynamical states of the Kuramoto model. 

As already mentioned, the four-dimensional system ( fTSl l is 
effectively three-dimensional due to the existence of a neu- 
tral global rotation. Interestingly the attractors of the model 
are apparently embedded into a fvvo-dimensional plane. Nu- 
merical simulations of Eqs. JTSl l using arbitrary initial con- 
ditions show that the dynamics always collapses into a plane 



which, by virtue of the neutral rotation e"^^, can be made co- 
incident with the (a;i, X2) plane, hereafter refeiTed to as the 
"real plane". The stability against perturbations transversal to 
the real plane (and not tangent to the global rotation) is dif- 
ficult to prove analytically. For the fixed point X(^2) born at 
the saddle-node bifurcation, the stability against transversal 
perturbations is proven in Appendix A. Other attractors (limit 
cycle) are transversally stable according to our numerical sim- 
ulations. 

Numerical simulations of the reduced Eqs. ( fTSl l with either 
real or complex variables, it is irrelevant, reveal that 

(i) The Hopf bifurcation at K = Kh is supercritical and 
it gives rise to a limit cycle around the origin. Due 
to the reflection symmetry of the equations z{t) van- 
ishes twice per period [this occurs when ai — ^a2, 
see Eq. (fTTbl. It is therefore reasonable to assume that 
the limit cycle coiTesponds to the SW state, for which 
the two counter-rotating clusters of phase-locked oscil- 
lators are tt out of phase twice per period. 

(ii) The oscillatory dynamics appearing at Kh is destroyed 
at i^T = Ksn where twin saddle-node bifurcations give 
rise to twin pairs of fixed points on the limit cycle. 
This bifurcation is known as SNIC (saddle-node on the 
invariant circle), or SNIPER (saddle-node infinite pe- 
riod). As K approaches Ksn from below the period 
of \z{t)\ diverges due to the slowing down of the dy- 
namics at the twin bottlenecks anticipating the cease of 
oscillations via the (double) SNIC bifurcation. 

Finally, numerical simulations of the full Kuramoto 
model ([T]l confirm the scenario I S W PS predicted by the 
reduced equations (fTsT l. We have numerically determined the 
boundaries of different behaviors: Square symbols in Fig. [3] 
are points in which the incoherent state loses stability leading 
to a SW state. Additionally, triangles indicate points where 
the order parameter becomes stationary. 



D. Concluding remarks 

Distribution (|5]l with ^ = 7 becomes unimodal only for 7 = 
0. As 7 — > the bimodal distribution tends to a unimodal, but 
the limit is nonregular. The remarkable point is that bistability 
is not observed, even if the central dip is extremely narrow 
(7^0). This is in sharp contrast with the scenario found 
when the peaks are close to merge with ^"(O) 0+ at the 
usual unimodal-bimodal transition (see below). 

Another interesting fact is that the counter-rotating clusters 
of the S W are born at the Hopf bifurcation ( l24b with frequen- 
cies although the maxima of the distribution are located 
at ±7. This means that the relative shift between distribu- 
tion's maxima and cluster frequencies at the onset of the SW 
diverges as 7 — > 0. This is a consequence of the extreme 
asymmetry of the peaks in this limit. 
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FIG. 4: Phase diagram for ^ — 0.5. Solid lines mark the bifurca- 
tions: Saddle-node off the limit cycle (SN), SNIC, Hopf bifurcation 
[Eq. l l24b 1. heteroclinic bifurcation (found numerically using the re- 
duced equations), and pitchfork bifurcation [Eq. | |30H . Three big 
circles signal the codimension-two points: (A) Takens-Bogdanov, 
(B) degenerate pitchfork, (D) saddle-node separatrix-loop. The open 
symbols correspond to different bifurcations found by numerical in- 
tegration of Eqs. O with — 2000. Filled symbols inside each 
region indicate parameter values for the phase portraits in Fig.[5l 



IV. BIMODAL DISTRIBUTIONS NONVANISHING AT 
THEIR CENTER < 7) 

In this section we analyze the case ^ < 7, which is com- 
plementary to the one studied in the previous setion = 7). 
Thus, in the present case we let ^ and 7 to be independent of 
each other (see Fig.|2]). As we did in the previous section, we 
determine first the local bifurcations of the fixed points, and 
then we summarize our findings in the (7, K) phase plane to- 
gether with the results obtained by numerical integration of 
the reduced Eqs. (fTSb as well as of the full Kuramoto model 



A. Fixed points 

1. The incoherent state and its stability 
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FIG. 5: Phase portraits in (rotated) xi,X2 coordinates for qualita- 
tively different cases. Each panel corresponds to a value of 7 and K 
at the position of a filled symbol in Fig.|4] (a,b) Partial synchroniza- 
tion with K = A, and (a) 7 = 0.6 and (b) 7 = 0.75; (c) Coexistence 
SW/PS: 7 = 0.67, K = 3.45; (d) Coexistence I/PS: 7 = 0.7, 
A' = 3.3 ; (e) SW, 7 = 0.6, it" = 3.5 ; (f) I, 7 = 0.6, K = 2.5. 



of Hopf and pitchfork bifurcations collide at the codimension- 
two point where Kh = Kp and ^ ~ ^a- This point is of the 
double zero eigenvalue type (Takens-Bogdanov) jsoll . 

The boundaries jTM and ( [30l l for Hopf and pitchfork insta- 
bilities have been also obtained following a different approach 
in Appendix B. 



2. Non-trivial fixed points (partial synchronization) 



The incoherent state becomes unstable in two possible ways 
depending on the value of with respect to: 



(29) 



(see line A in Fig.|2|i. For ^ < £^a, there is a degenerate Hopf 
bifurcation at the critical value Kh given by Eq. ( |24] | which is 
independent of ^. For ^ > £^a, the instability of the incoherent 
state occurs via a pitchfork bifurcation at: 



Kp = 



7-e ■ 



(30) 



The bifurcation is subcritical, and it switches to supercritical 

1 /3 

when the distribution becomes unimodal at 7 > . The loci 



A saddle-node bifurcation occurs when P(X) in Eq. ( |20l ) 
has exactly two roots (one of them two-fold degenerate). And 
this bifurcation point can be determined numerically finding 
the value of k where the discriminant of P{X) vanishes. The 
scenario is similar to the one observed for ^ = 7, but in this 

(1) > exists up to the pitchfork 
Kp. If the distribution is 



case the saddle solution X, 



bifurcation with the origin at K 
unimodal < what makes this solution not valid. 



B. Numerical simulations and phase diagram 

Our analytical results provide information about local bifur- 
cations. In addition we have performed numerical simulations 
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of the ODEs ( fTsT l. in order to obtain the full system's picture. 
As occurred in the previous section, we can assume that OLj 
are real variables. In addition, we have performed numeri- 
cal simulations of the original system that indicate that this 
assumption yields to correct results. 

Figure|4]shows the disposition of qualitatively different dy- 
namics in the parameters space spanned by 7 and K, for a par- 
ticular value of ^. Like in 1 13] we find that three codimension- 
two points organize the parameter space: Takens-Bogdanov 
(A), degenerate pitchfork (B), and Saddle-node separatrix- 
loop (D) llsill . The three codimension-two points collapse at 
^ = 7 = 0, see Fig. |2] and expressions ( |29] l and (|7]l. Line 
D approaches the origin linearly: ^d(7 ^ 0) = a7 with 
a ~ 0.493, suspiciously close to i. 

One can better understand Fig. |4] looking at the panels 
of Fig. |5] in which phase portraits for qualitatively different 
states are shown. In the rightmost part of Fig. S] 7 > 75 = 
the distribution becomes unimodal, and thus the stan- 
dard route to partial synchronization is found. In the leftmost 
part, C < 7 < 7D - 0.59997 {Kd ^ 3.7646), we have 
the same route than in the previous section, i.e. a SW state 
limited by Hopf and SNIC bifurcations. In contrast, in the 
central part of the phase diagram (around point A), there ex- 
ist two regions with bistability where the observed asymptotic 
state depends on the initial conditions. In one region (SW/PS) 
standing waves and partial synchronization coexist, and the 
SW state (a limit cycle) disappears via a heteroclinic collision 
with the saddle points born at mirror saddle-node bifurcations. 
In the second region (I/PS) incoherence and partial synchro- 
nization coexist. 

Bifurcation lines in Fig.|4]are calculated from analytical re- 
sults and from numerical integration of the ODEs ( fTsT l. Empty 
symbols in the figure show the bifurcations determined inte- 
grating the Kuramoto model with N = 2000. The agreement 
is good and confirms the validity of the OA ansatz. 



1. Codimension-two point A 

In this subsection we make a short digression about the 
codimension-two point A and the importance of the symme- 
tries in the model. Point A in Fig. |4] is a Takens-Bogdanov 
point of system (fTsT l that has 0{2) symmetry. This stems 
from the inherent 0(2) symmetry of the Kuramoto model 
[with symmetric g{uj)]. Numerics show that the asymptotic 
dynamics occurs in the real plane — i.e. Eqs. (fTSl l with real 
coordinates — where the symmetry group is only Z2 C 0(2). 
This symmetry imposes the global bifurcation (Het) to be non- 
tangent to the Hopf line |30], in contrast with a nonsymmetric 
Takens-Bogdanov point. Two scenarios are possible around 
the odd-symmetric Takens-Bogdanov point IBOll . Hence, one 
may wonder if the alternative scenario, involving a saddle- 
node bifurcation of limit-cycles, might also be found in the 
Kuramoto model. 

The scenario that we have presented in this section (see also 
ifisll ) is apparently the same one Bonilla et al. lITsIl uncovered 
in the neighborhood of the Takens-Bogdanov point for the Ku- 
ramoto model with additive noise and a bi-delta frequency dis- 



tribution. In that work the full 0(2) symmetry is taken into 
account. Refs. [12., 15] found that, due to the 0(2) symmetry, 
the degenerate Hopf bifurcation gives rise to a branch of un- 
stable traveling wave solutions, in addition to the stable SW. 
According to [15J these traveling wave solutions should dis- 
appear at a certain K < Kp in a local bifurcation with the 
saddle fixed points born at the SN bifurcations. This bi- 
furcation reverses the transversal stability of the saddle fixed 
points, what in turn makes congruent the pitchfork bifurcation 
of these fixed points with the completely unstable fixed point 
at origin. We think these traveling wave solutions and their 
associated bifurcations are captured by the reduced Eqs. ( fTSb 
because the OA ansatz has retained the 0(2) symmetry of the 
model. This means that although the relevant dynamics (the 
attractors) are inside the real plane of (ai, 02), physical un- 
stable objects (traveling waves) "live" outside this plane. 



V. CONCLUSIONS 

We have investigated the routes to synchronization in the 
Kuramoto model with a bimodal distribution constructed 
as the difference of two unimodal distributions of different 
widths. These distributions admit an arbitrarily deep and nar- 
row central dip, what is not achievable in distribution types 
considered in the past. This has allowed us to reinforce and 
extend the results recently published in [13]. 

We have found that bimodal distributions (|5]) near uni- 
modahty produce hysteretic phase transitions, except in some 
region in the neighborhood of the unimodal limit (^, 7) = 
(0,0), see Fig. 12 

We expect a wide family of bimodal distributions to ex- 
hibit the same qualitative features that Fig. |2l The hysteretic 
region exist at the bimodal side of the unimodal-bimodal bor- 
der, and it shrinks as the nonregular unimodal-bimodal tran- 
sition (^"(0) = 00) is approached. Moreover the absence of 
hysteresis for g{Q) = should be found in any bimodal distri- 
bution if the dependence is quadratic — as in our distribution 
(|5]l — or has a larger power: g{Ld) cx \loY' for small cj, with 
1/ > 2. 
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APPENDIX A: PROOF OF THE TRANSVERSAL 
STABILITY OF FIXED POINT X(2) IN SEClml 

Global phase shift invariance, (ai,Q;2) (aie'^, a2e''^), 
allows to reduce Eqs. ( fTSl l in one dimension by passing to 
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polar coordinates, aj = pje^'^^ , and defining the phase differ- 
ence tp = (f>i — (f>2- We obtain three ODEs: 



Pi = -pi + Hpi - ^P2Cosi/))(l - pI) 



P2 = -IP2 + k{Pi costp - CP2)(1 - pI) 



PiP2^ = -k [(1 - ^)pIpI + pI 



ipl] sin V 



(Ala) 



(Alb) 



(Ale) 



In Sec. ini] we took ^ = 7 and found that twin saddle- 
node bifurcations (namely SNICs) give rise to two pairs of 
fixed points. Here we prove (we rather sketch the proof) the 
transversal stability of the mirror fixed points associated to 



X, 



lEq. (EB. 



(2) 

First of all note that ^(2) yields a fixed point (xi, 2:2), and 
its mirror image, with xi and X2 having the same sign, -0 — 0. 
This is a consequence of Eq. ( |2T] ) because X(2) < 1 — 1/fc. 
The latter inequality stems from the fact that Q{1 — 1/fc) = 
—7^ < and by continuation of the solutions from fc = 00: 
limfe^oo -'^(i) (fc) = 0, limfe^oo -'^(2,3) {k) = 1. 



Therefore we have to prove that factor 



F={l--f)pipi 



pI 



ipI 



(A2) 



in Eq. ( I A let for ij} is positive. Replacing p\ = X^2) p\ 
X(2) + 1/fc, we have 



F = (1 - 7) [Xf^) + X(2) (1 + 1/fc)] - 7/fc- 



(A3) 



As X(2) exists only above the saddle-node bifurcation (fc > 

ksN) and ksN > fc/f = (1 + 7)/(l - 7). 



with 



F > (1 -7)/i 



(A4) 



(A5) 



Then > is a sufficient condition for the transversal stabil- 
ity of the fixed point. 

It suffices to prove that h is positive at the locus of the 
saddle-node bifurcation because X(2)(fc;7) exhibits its min- 
imal value over fc precisely at the bifurcation: X(2)(fc > 
ksN, 7) > X(2) (fcsjv, 7)- For our aim it is better to parame- 
terize the SNIC line by k instead of 7. Hence we to introduce 
in iA5i the expressions 



(i) 7 as a function of kg^, via Eq. 



(ii) X(^2){ksN), determined from dZST l in the two-fold root 
case. 

The calculation of terms (i) and (ii) can be readily done with 
symbolic software such as MATHEMATICA. As a result we 
obtain a function h{ksN) that is positive in all the domain of 

ksN e (1,00). 

Moreover using expressions (IZTT i and (|28T l we can get ap- 
proximate expression for h (as a function of 7): 



M7-0)= (|)'^' + 0(7) 



(A6) 



^ ^ ^ 2/3 

h{j -> 1) ~ 0.0858 (A7) 

APPENDIX B: STABILITY OF THE INCOHERENT STATE 
IN THE KURAMOTO MODEL WITH NOISE 

For the sake of completeness, and as a double-check of 
some of the results obtained, we study here the stability of 
the incoherent state when the model is perturbed with addi- 
tive white noises. In this case, the right hand side of Eq. ([T]) 
has to be supplemented with uncorrected fluctuating terms rji 
satisfying {'rii{t)rij{t')) — 2a5ij5{t ~ t'). So far a counterpart 
of the Ott-Antonsen ansatz for the stochastic problem has not 
been found. It is nonetheless possible to obtain the stability 
boundary of incoherence resorting to the Strogatz and Mirollo 
relation for the discrete spectrum of eigenvalues A 12 ill : 



K 



2 X + (T + ioJ 



-duj = 1. 



(Bl) 



Considering the distribution of frequencies (|5]), this equation 
can be solved for the eigenvalues A. 

Noise increases the domain of the incoherent state. Hopf 
and pitchfork bifurcations continue to occur, but the values of 
K are shifted to larger values. We obtain: 



Kh = 2 + 27 + 4(7 

2(7 + a)(l-0(l + ^) 



Kp = 



(7-0+^(1-0 



(B2) 
(B3) 



that indeed reduce to Eqs. (l24l i and dSOl l for a = 0. The loca- 
tion of the Takens-Bogdanov point [c.f. Eq. ( |29l )l also varies 
and now pitchfork and Hopf bifurcations collide (Kh = Kp) 
at : 



u = 



1 - 



(B4) 
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